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We consider the hydrodynamic modes for dense trapped ultracold gases, where the interparticle 
distance is comparable to the scattering length. We show that the experimental determination of 
the hydrodynamic mode frequencies allows to obtain quite directly the equation of state of a dense 
gas. As an example we investigate the case of two equal fermionic populations in different hyperfine 
states with attractive interaction and in particular the vicinity of the collapse. 



Most of the fascinating recent work on ultracold gases |^] has been dealing with dilute situations. Naturally even 
in this regime interactions play an important role, as in the case of BEC where they strongly increase the size of 
the condensate compared to the free boson case. In this dilute regime the scattering length is small compared to 
the interparticle distance. However it is of great interest to explore the dense gas regime where scattering length 
and interparticle distance are comparable. This would lead to physical systems which are very simple examples of 
strongly interacting systems. These have much more complicated counterparts in condensed matter physics, such as 
liquid ''He or ■^He, or the electron gas in metals. This regime is also of major experimental interest in the search for a 
BCS superfluid in fermion gases ^J, since this is the range where the higher critical temperatures will be found, 
which should make the transition more accessible. Experimentally this dense regime corresponds to large scattering 
lengths, which can be reached in the vicinity of Feshbach resonances, as it has already been seen in optical traps j^. 

In this paper we show that the experimental determination of the mode frequencies in the hydrodynamic regime 
allows to obtain quite efficiently and directly the equation of state of a dense gas. Hydrodynamic equations are valid 
in the limit of low frequency and long wavelength, and they have already been used to study the dilute Bose gas 
(with very good agreement with experiment) as well as the free Fermi gas Q . In these cases the equation of state is 
known. We show that, rather surprisingly, the analysis of the equation giving the mode frequencies is not much more 
complicated when the equation of state is unknown and that one can conveniently invert the problem and obtain 
the equation of state from the mode spectrum. Although one might consider the extension to higher temperature 
we work in the low temperature range where thermal effects are small and we neglect dissipation, so we deal with 
a perfect fluid. This should be valid for a strongly degenerate Fermi gas in its normal state. Naturally our results 
apply also to a superfluid when normal liquid effects are negligible, such as low temperature Bose condensates or low 
temperature simple scalar BCS superfluids. We consider for simplicity an isotropic trapping potential V{r) (mostly 
the harmonic case) but, somewhat surprinsingly most of our procedures can be generalized to the case of anisotropic 
harmonic traps. Also we treat the 3D case, but lower dimensions can be handled in the same way. As an example 
we apply our treatment to the case of two equal fermionic populations in different hyperfine states with attractive 
interaction, in particular we investigate the vicinity of the collapse, a very interesting physical situation analogous to 
■^Li BEC collapse. 

With our hypotheses hydrodynamics reduces to Euler equation mndv/dt = —VP — n\7V{r) supplemented by 
particle conservation dn/dt + V{nv) = for density n(r,t). The equilibrium particle density no(r) satisfies ^(no(r)) + 
V{r) — jl where /2 is the overall chemical potential. Below we refer for short to ^{n) as the equation of state. Linearizing 
these equations around equilibrium one finds that the density fluctuation ni(r, i) = n(r,t) — no(r) oscillating at 
frequency uj satisfies V^(nii9P/9no) + V(niVy) + muu'^ni = with dP/dno = nod^i/dnQ. Of particular interest 
is the 'neutral mode' solution ^^(r), corresponding to the density fluctuation produced by a small shift Sp, of the 
overall chemical potential, that is n5(r) = {diJ,/dno)^^dfi. Since the result is still an equilibrium situation, this mode 
corresponds to a; = (as can be checked directly). Actually this mode does not correspond to a physical mode since 
it does not conserve particle number, whereas this number is automatically conserved by all the solutions with non 
zero frequencies, as it can be seen by integrating the starting equations over the whole gas volume. 

We make the change ni(r) = ni{r)w{r)Yim (i-e. the local fluctuation of the chemical potential is a convenient 
variable). Since w^r) = 1 corresponds to the neutral mode, we may expect a simple equation. Indeed we obtain : 



PACS numbers : 03.75.Fi, 32.80.Pj, 47.35. +i, 67.40.Hf 



rw" + [2 + rL\r)]w' - [ 



L'{r)\w 



(1) 



+ 



V'{r) 



r 



1 



where we have set L{r) = ln(no(r)) with L'{r) = dL/dr, and V'{r) = dV{r)/dr. This equation for the mode 
frequencies has a quite simple form. In particular, as soon as V{r) is known, the properties of the fluid appear 
only through the function L{r), which is itself simply related to the equation of state through /i(no(r)) + V{r) = fl. 
Therefore it appears much more convenient to model L{r), rather than Indeed Eq.(l) lends itself to a very 

large number of specific models with analytical solutions or quasi-analytical solutions, as we will see below in the 
case of the harmonic trap. Before proceeding to this case it is also interesting to note that Eq.(l) may be written 
with the form of a one-dimensional Schrodinger equation (with energy equal to zero) by making the further change 
w{r) = ip{r)l{r^JnQ{r)). The corresponding potential is found to be {\ / r^muj'^ /V')L' +{\ / A)L''^ + {1 /2)L" + /r'^ 
and is again simply related to L{r). This form is of particular interest when one has an explicit analytical solution 
for an approximate model, as we will find below. One can then easily correct the results by a standard first order 
perturbation calculation. However Eq.(l) is more convenient for actual solutions. 

Let us now specialize to the case of the harmonic trap V(r) = ^mil^r^. It is then convenient to make the further 
change w{r) = r^v{r) to get rid of the l{l + \)/r term. This leads to : 

rv" + [2{l + 1) + rL'{r)]v' - {v^ - l)L'{r)v = (2) 

where = up' . We check on this equation that, whatever the equation of state of the fluid, we have as 
expected the dipole mode (/ = 1) solution corresponding to the oscillation in the harmonic trap of the gas as a whole, 
at frequency = fJ. It corresponds |^ to the solution v — \ with = \. Furthermore v — \ gives also uo = r2\/I 
whatever L{r) i.e. independent of the equation of state of the fluid (and in particular whether it is a Bose or Fermi 
gas) This generalizes for an interacting fluid, at low temperature, results obtained by Griffin et al. and Stringari 
for a Bose gas and by Bruun and Clark ||^ for a free Fermi gas. Naturally one can check on the starting equations 
that the corresponding 'surface mode' |9| fluctuation n\(r) ~ r^rb{(r) cx r''~'^ [driQ / dr) is indeed solution. 

Let us note first that Eq.(2) is invariant under the change of scale r — > Kr, provided we make the same change of 
scale for L{r). This allows us to take in the following the gas radius R as unity (consistently with hydrodynamics we 
work within the Thomas- Fermi approximation). Next we notice that Eq.(2) is only slightly modified by the change 



of variable y = r" provided again that the same change is made for L{r). This gives |11 
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A convenient feature of Eq.(3) is also that the absolute scale in density disappears in L' and only n{r) = n{r)/n{0) 
enters. We introduce similarly a normalized local chemical potential fl{r) = fi{n{r)) / fi(n{0)) where /i(n(0)) is simply 
obtained from the gas radius R by /i(n(0)) = ^m^l^R^, leading to ft — 1 ~ P, with r being now in units of R. 

Looking now for sirnple situations where we can solve Eq.(3), we consider naturally first the case of the non 
interacting Fermi gas |^| . This gives p, = n^^^ with p = 3/2. Similarly we can consider an interacting dilute Bose 
gas Q where ji — gn {g being the coupling constant) leading again to [1 = n^/P with now p ~ 1. These two cases 
imply L{r) = pln(l — r")with a = 2 in Eq.(3). In this way we are lead to consider for any a and p the model 
dL/dy = — p/(l — y) for which Eq.(3) becomes explicitely : 

yi^^v)-r2+[c~v{P + c)]— +p v^Q (4) 

dy^ dy a 

with c = 1 + The general solution of this equation, giving a non divergent density fluctuation for r = 0, is 

the hypergeometric function F{a, b; c; y), with a + b= p + c— 1 and ah = —p '^~^ . We have furthermore to require 
that the solutions satisfy the boundary condition that the outgoing particle current is zero everywhere on the sphere 
r = 1. This is not verified |p!^ by the general solution, except if we require that a = —n where n is a non negative 
integer in which case the solution is a polynomial Q . This leads to the following general result for the normal mode 
frequencies : 

, a , 21 + 1. 
= l + -n{n+p+ ) (5) 

which agrees naturally with Stringari Q for a = 2 and p — I, and with Bruun and Clark for a — 2 and p — 3/2. 
One may naturally wonder about the interest of these results for other values of our parameters a and p. These cases 
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correspond to the density distribution n(r) = (1 — r")^ and the equation of state /2 = 1 — (1 — n^/^*)^/". Our point is 
that these corresponding models can be used to represent closely the equation of state fi{n) for a general fluid (with 
a given maximum density n(0)). We will show explicitely below that the flexibility offered by the two parameters a 
and p makes it a very convenient and efficient procedure. However at first glance these general models do not seem 
very physical since, although their density is properly vanishing at the gas radius, they give near this border p, « n^^^ 
whereas one should get the dilute gas behaviour p = 1 (bosons) or p = 3/2 (fermions). However just because the 
gas is dilute near r = 1, we do not expect this part of the gas to play a significant role. Similarly these models give 
n « 1 — pr" for small r whereas one expects only the case a = 2 to occur for a regular equation of state. But if 
despite of these shortcomings no{r) is closely approximated over the whole range, one may expect this modeling to 
be already quite reasonable. 

Nevertheless it is clearly of interest to consider more complicated models which could display proper behaviour 
near the center and the border of the cloud. Although we have not obtained such models with completely analytical 
solutions, we have found a large class of models with quasi-analytical solutions which are in practice not different 
from fully analytical solutions. These are the models dL/dy — — J2k=oP'^y^ /(^ ~ v) (where we could take a = 2 and 
p = X]fc=o — 1 or 3/2 in order to have the proper center and border behaviour). To be simple and specific let us take 
the case K = 1, giving — dL/dy = [po +pxy)/{l — y). This corresponds to the equation of state n = jF exp[pi(l — Jl)]. 
In this case Eq.(3) becomes : 



2/(1 - y)^ + (g2y +qiy + qo)-^ + {riy + ro}v = (6) 



with ~ ^Pi, Qi = ^(1 + +Po)i 90 = 1 + ''I = Pi^^-—^ and tq — po^^—^- When we look for a series expansion 
of the solution v = "^n^QCLny", we find the following recursion relation (with a_i = 0) : [{n + l){n + (7o)]'^n+i + 
[—n(n — 1) + Jiqi + ro]a„ + [{n — l)q2 + ri]an-i = which does not allow in general for a polynomial solution. For large 
n this relation becomes asymptotically an+i — a„ = —{q2/n)an-i- This leads to the standard behaviour a„+i w an 
giving a convergence radius equal to 1 . This is the same situation as for the hypergeometric function in Eq. (4) (which 
corresponds to q2 = ri = 0) and this leads in the same way to a singular behaviour for y = 1 which does not agree 
with the boundary condition. On the other hand the above asymptotic relation may also have solutions Un+i <C a„ 
implying a„ « ('Z2/^)«n-i which gives a„ ~ 1/n!. This very rapidly convergent series has an infinite convergence 
radius and no singularity for j/ = 1. It corresponds to the physically acceptable solutions. Since we have only to 
deal with the range y € [0, 1] this solution is a quasi-polynomial since the higher order terms in the series are very 
rapidly negligible. This is quite analogous to the polynomial solution of the hypergeometric differential equation. 
Naturally these solutions arise only for special values of our parameters, which gives finally the mode frequencies. 
In practice these parameters are found very easily in the following way. We solve iteratively the recursion relation 
for a„ with < n < A'^ and we require a^r+i = (as if we had a polynomial solution). Since vq and ri are linear 
in i/^ — I, this is equivalent to find the roots of an equation of order N for ~ I. We then increase the value of 
N = 1, 2, ... until the roots have converged. For the lowest root this is usually a very fast convergence, so one could 
obtain approximate analytical expressions. But the numerics is so easy that this seems unnecessary. Naturally this 
procedure gives immediately Eq.(5) for the hypergeometric case. All this analysis and procedure can be extended to 
the case of A" > 1, which gives a„ ~ 

(l/n!)i/^ for the convergence of the series for the physical solution. 

As an example we turn now to a specific case and consider the case of two equal populations of fermions in different 
hyperfine states. This may be the case of ^Li or "^^K near a Feshbach resonance |10|. We assume an attractive 



interaction between unlike atoms with an interaction related to the (negative) diffusion length hy g — 47r?i^a/m, 
and we take the Hartree approximation to describe this system. For total atomic density n, the chemical potential is 
given by ^(?^) = h^kp/2m — \g\n/2 with 37r^n = kp. More generally we could deal with situations where /i(n) is a 
polynomial in u = kp. In an isotropic harmonic trap the equilibrium density satisfies then = y ~ P{u) where P{u) 
is a polynomial. In order to solve directly this case it is more convenient to rewrite Eq.(2) (taking a — i) with the 
variable u. This leads to: 

PP'v"+ [{l + l)p'^+-PP'^pp"]v'- 3(^'-0 :Pl^^ 
2 u 2 u 

with P' = dP/du and P" — d^P/du^. Since this equation is homogeneous in u we can rescale this variable and have 
it varying between and 1. Similarly the homogeneity in P allows to conveniently rescale it. The free fermion case 
corresponds to P{u) = 1 — m^. After taking as a variable one obtains the hypergeometric differential equation 
which leads again to the free fermion result m |(i/^ ^ = "-("- + ^ + 2). For the Hartree approximation we can write 
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P{u) = 1 — — §-^(1 ^ '^"^) with the couphng constant A = 2fci?(0)|a|/7r and A:f(0) the equihbrium Fermi wavevector 
at the center r = 0. This couphng constant goes from 0, for the very dilute regime, to 1 when we reach at the center 
the instabihty where the gas is going to coUapse under the overaU attractive atomic interaction. In this case we can 
not solve Eq.(7) by quasi- polynomials as Eq.(6) because this equation has five regular singular points, instead of only 
two for Eq.(6). 

We have solved Eq.(7) numerically, as a function of A, for the first three monopole mode frequencies {1 = 0). The 
results are given in Fig.l . One sees that the frequencies decrease for increasing attractive interaction. This was to 
be expected since the gas gets more compressible when one approaches the instability. However, in contrast with 
what one might anticipate, we do not find the lowest mode frequency going to zero at the instability. This can 
be understood because the instability density is reached only at the center, and the rest of the gas still provides a 
restoring force accounting for the nonzero frequency. Actually this limit can not be reached experimentally since 
the modes correspond to an infinitesimal density oscillation. Any finite oscillation will induce nonlinear effects and 
produce a collapse of the gas. So one should experimentally work with ever smaller oscillation when one goes near 
the instability. We believe also that including quantum effects (the hydrodynamic description is not correct at the 
scale of the Fermi wavelength) will, so to speak, smear the region where the instability density is reached and lead to 
a zero frequency mode very near the instability (preliminary calculations support this view). 

It is now of interest to consider an approximate solution of this same problem with the modeling we have discussed 
above. With the Hartree approximation the relation between the reduced chemical potential jl and the reduced 
density n is /x = (3n^/^ — 2An)/(3 — 2A). For each value of A we approximate p,{n) by /i = 1 — (1 — n^^^)^/", where 
we obtain the parameters p and a by a least square fit. Then the mode frequencies are given by v"^ = ^n{n + p + ^) 
for n = 1,2,3. The results are given in Fig.l and they are surprisingly close to om* exact results from numerical 
integration, although our modeling does not give in general a good description of fi{n) neither for n nor for n ~ 1, 
as already mentioned. The real interest of this approximate treatment is that it is very easily inverted and allows to 
analyze readily experimental data. If we have, for a given A, the frequencies of the first two modes, we can obtain the 
values of p and a from the formula for the frequencies (the values for the higher modes give a check on this modeling). 
Hence we have p.{n) from our model. By varying the coupling constant we can obtain the equation of state /x(n) 
it originates from. In our case we could, from the result of Fig. 1 , recover the Hartree expression for /x(n) to a very 
good approximation. This method is much more flexible than trying to fit the data with a specific equation of state. 
Naturally it is not clear that this modeling will always work so well and we can make use of the richer models with 
quasi-polynomial solutions we have already discussed, as we will see below. 
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FIG. 1. Reduced mode frequency = [lu/Q)^ for a Fermi gas within the Hartree approximation as a function of the coupling 
constant A. Full line: exact result from the numerical solution of Eq.(7). Dashed line: approximate analytical solution. Insert: 
lowest mode frequency at the collapse (A = 1) for a generalized equation of state (see text). The straight line is = 1 + 2P„. 
Pk = 0.577 for the Hartree approximation. 

In this aim we investigate the situation when the gas is right at the collapse limit. This corresponds to A = 1 in 
the Hartree approximation. However it is not clear that the Hartree approximation will be so good, in particular in 
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the vicinity of the collapse and we may consider more general equations of state. The critical density at the collapse 
Uc satisfies d^/ dn^ = 0. Let us assume that ^{n) is regular in the vicinity of Uc, so we can write the expansion 
— iJ,c — n{n — nc)^ /2. At the collapse we have n{r = 0) = ric and /ic — ■^mil^R^ where Rc is the gas radius 
at the collapse. In the vicinity of the center + ^mil^r^ — fic gives 1 — n{r) — {ric — n)/nc — {r / Rc){Rc/ LfS)^ 
where we have introduced the length scale = {nc/fl){K/my/^ . This behaviour can be described by our model 
n(r) = (1 — r")P (with now r in units of Rc ) if we take a — 1 and p = Rc/Li^. Then from Eq.(5) the lowest 
frequency monopole mode is given hy = 1 + 2/p. The Hartree approximation corresponds to p = \/3- The result 
1/^ = 1 + 2/-\/3 = 2.15 is in excellent agreement with Fig.l (actually our least square fit gives a and p very near 1 
and 73). If we had chosen p — 3/2 in order to have a better description of near r — 1, the result i^^ = 7/3 would 
not be as good. Now we can improve our modeling by taking h(r) = (1 — r)'^/^ exp(pir) which has quasi-polynomial 
solutions. We have already taken our two parameters a = 1 and p = 3/2 to have a good description near the center 
and the border. We can still adjust pi to have a better description n{r) — 1 ~ rRc/L^ near the center. This gives 
Pi = 3/2 — RcjLf^. Our results for the lowest monopole frequency are reported in the insert of Fig.l as a function 
of the parameter = Li^/Rc = (Kn^/2/Zc)^/^. The results of the two models are almost indistinguishable in the 
vicinity of the Hartree case (P„ = l/-\/3) and go to the same limit — 1 for (see below). In between they 

differ somewhat. It is clear that n(r) = (1 — rY is not flexible enough to give a precise description of n(r), whereas 
n(r) = (1 — r)'^/^ exp(pir) satisfies all the constraints we have set, so we can be fairly confident in the results. We 
see in Fig.l that the experimental determination of the lowest monopole frequency gives direct information on the 
equation of state near critical density since it gives immediately P^ — [Kric/'^iicY^'^ . Let us also mention that the 
monopole modes fluctuations are divergent for r = in this collapse limit, as one could expect physically. 

The limit where the parameter P^ goes to zero is particularly interesting. It corresponds to the situation where 
most of the gas is concentrated at small r and the rest of the gas is irrelevant. One may take Rc ^ oo and take 
Lk as new length unit for r. Then our model n(r) = (1 — r)'^/^ exp(pir) becomes n(r) = exp(— r), corresponding to 
p, = 1 — In^(n). The general result is then easily obtained from Eq.(2) : since we have merely L'{r) = —1, Eq.(2) 
reduces to the equation for Laguerre polynomials, and the spectrum for all the modes is degenerate and given by 
v"^ = N where A'^ is a positive integer. Actually the corresponding Schrodinger form of the equation, mentionned 
above, reduces to the one for the hydrogen atom. Similarly in our model n[r) = (1 — r")P, we have p — > cxd in this 
limit. When we take as new length unit for r, this model becomes again n{r) — exp(~r) with the same result for 
the spectrum. Indeed we check on Eq.(5) that, in the limit p ^ oo, the spectrum is given hy = I + n = N and is 
degenerate. 

We are very grateful to M. Brachet, Y. Castin, C. Cohen- Tannoudji, J. Dalibard and C. Salomon for very stimulating 
discussions. 
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